Biostat/Biomath M257
This lecture introduces GPU computing in Julia.
GPUs are ubiquitous in modern computers. Following are NVIDIA GPUs on today’s typical computer systems.
| NVIDIA GPUs | H100 PCIe | RTX 6000 | RTX 5000 |
|---|---|---|---|
 |
 |
 |
|
| Computers | servers, cluster | desktop | laptop |
 |
 |
 |
|
| Main usage | scientific computing | daily work, gaming | daily work |
| Memory | 80 GB | 48 GB | 16 GB |
| Memory bandwidth | 2 TB/sec | 960 GB/sec | 576 GB/sec |
| Number of cores | ??? | ??? | ??? |
| Processor clock | ??? GHz | ??? GHz | ??? GHz |
| Peak DP performance | 26 TFLOPS | ??? TFLOPS | ??? TFLOPS |
| Peak SP performance | 51 TFLOPS | 91.1 TFLOPS | 42.6 TFLOPS |
GPUs contain 1000s of processing cores on a single card; several cards can fit in a desktop PC
Each core carries out the same operations in parallel on different input data – single program, multiple data (SPMD) paradigm
Extremely high arithmetic intensity if one can transfer the data onto and results off of the processors quickly
 |
 |
|---|---|
 |
 |
GPU support by Julia is under active development. Check JuliaGPU for currently available packages.
There are multiple paradigms to program GPU in Julia, depending on the specific hardware.
CUDA is an ecosystem exclusively for Nvidia GPUs. There are extensive CUDA libraries for scientific computing: CuBLAS, CuRAND, CuSparse, CuSolve, CuDNN, …
The CUDA.jl package allows defining arrays on Nvidia GPUs and overloads many common operations.
The AMDGPU.jl package allows defining arrays on AMD GPUs and overloads many common operations.
The Metal.jl package allows defining arrays on Apple Silicon GPU and overloads many common operations.
AppleAccelerate.jl wraps the macOS Accelerate framework, which provides high-performance libraries for linear algebra, signal processing, and image processing on Apple Silicon CPU. This is analog of MKL for Intel CPU.
The oneAPI.jl package allows defining arrays on Intel GPUs and overloads many common operations.
I’ll illustrate using Metal.jl on my MacBook Pro running MacOS Sequoia 15.4. It has Apple M2 chip with 38 GPU cores.
versioninfo()Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 24 × 13th Gen Intel(R) Core(TM) i7-13700
WORD_SIZE: 64
LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 1 default, 0 interactive, 1 GC (on 24 virtual cores)Load packages:
using Pkg
Pkg.activate(pwd())
Pkg.instantiate()
Pkg.status() Activating project at `~/biostat-257-25S/slides/09-juliagpu`Status `~/biostat-257-25S/slides/09-juliagpu/Project.toml`
[6e4b80f9] BenchmarkTools v1.6.0
[bdcacae8] LoopVectorization v0.12.172
[8f75cd03] oneAPI v2.0.2
[37e2e46d] LinearAlgebra v1.11.0using oneAPI
oneAPI.versioninfo()Binary dependencies:
- NEO: 24.26.30049+0
- libigc: 1.0.17193+0
- gmmlib: 22.3.20+0
- SPIRV_LLVM_Translator_unified: 0.7.1+0
- SPIRV_Tools: 2025.1.0+1
Toolchain:
- Julia: 1.11.5
- LLVM: 16.0.6
1 driver:
- 00000000-0000-0000-1838-5f6401037561 (v1.3.30049, API v1.3.0)
1 device:
- Intel(R) Graphics [0xa780]using Random
Random.seed!(257)
# generate SP data on CPU
x = rand(Float32, 3, 3)
# transfer data form CPU to GPU
xd = oneArray(x)3×3 oneArray{Float32, 2, oneAPI.oneL0.DeviceBuffer}:
0.145793 0.939801 0.479926
0.567772 0.577251 0.81655
0.800538 0.38893 0.914135# generate array on GPU directly
# yd = Metal.ones(3, 3)
yd = oneArray(ones(Float32, 3, 3))3×3 oneArray{Float32, 2, oneAPI.oneL0.DeviceBuffer}:
1.0 1.0 1.0
1.0 1.0 1.0
1.0 1.0 1.0# collect data from GPU to CPU
x = collect(xd)3×3 Matrix{Float32}:
0.145793 0.939801 0.479926
0.567772 0.577251 0.81655
0.800538 0.38893 0.914135using BenchmarkTools, LinearAlgebra, Random
Random.seed!(257)
n = 2^14
# on CPU
x = rand(Float32, n, n)
y = rand(Float32, n, n)
z = zeros(Float32, n, n)
# on GPU
xd = oneArray(x)
yd = oneArray(y)
zd = oneArray(z);# SP matrix dot product on CPU: tr(X'Y)
bm_cpu = @benchmark dot($x, $y)BenchmarkTools.Trial: 88 samples with 1 evaluation per sample. Range (min … max): 54.340 ms … 59.727 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 57.093 ms ┊ GC (median): 0.00% Time (mean ± σ): 57.037 ms ± 1.232 ms ┊ GC (mean ± σ): 0.00% ± 0.00% ▂▂ ▂█ █ ▂ █▅ █ ▂ ▅▁█▁▁▁▅▁█▁▁▅▁▅███████▅▅██▅█▅▁███▅███▅██▁█▅█▁▅▁▅█▁▁███▅▅▁▁▁▅ ▁ 54.3 ms Histogram: frequency by time 59.6 ms < Memory estimate: 0 bytes, allocs estimate: 0.
# SP matrix dot product on GPU: tr(X'Y)
# why are there allocations?
bm_gpu = @benchmark oneAPI.@sync dot($xd, $yd)BenchmarkTools.Trial: 10 samples with 1 evaluation per sample. Range (min … max): 529.678 ms … 533.094 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 531.661 ms ┊ GC (median): 0.00% Time (mean ± σ): 531.526 ms ± 1.059 ms ┊ GC (mean ± σ): 0.00% ± 0.00% █ █ █ █ █ █ ██ █ █ █▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁█▁▁▁█▁▁▁▁▁▁█▁▁▁▁▁▁▁█▁▁▁██▁▁▁▁▁▁▁█▁▁▁▁▁▁▁█ ▁ 530 ms Histogram: frequency by time 533 ms < Memory estimate: 153.11 KiB, allocs estimate: 1877.
# speedup on GPU over CPU
median(bm_cpu.times) / median(bm_gpu.times)0.10738569883334086# SP broadcast on CPU: z .= x .* y
bm_cpu = @benchmark $z .= $x .* $yBenchmarkTools.Trial: 54 samples with 1 evaluation per sample. Range (min … max): 89.552 ms … 100.786 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 93.440 ms ┊ GC (median): 0.00% Time (mean ± σ): 93.341 ms ± 1.816 ms ┊ GC (mean ± σ): 0.00% ± 0.00% ▁ ▁▄ ▁ ▁ ▁▁▁ █ ▄ ▁ ▁ ▆▁▆▁▁▁▁▁▆▆▁▁▁▆▆▁▁▆▁█▁▁▆▁▆▆██▆█▁▆▆▆▆█▁███▆▆▁█▆▁█▆▁▆▆█▆▆▁▆▁▆▁█ ▁ 89.6 ms Histogram: frequency by time 95.8 ms < Memory estimate: 0 bytes, allocs estimate: 0.
# SP broadcast on GPU: z .= x .* y
# why is there allocation?
bm_gpu = @benchmark oneAPI.@sync $zd .= $xd .* $ydBenchmarkTools.Trial: 43 samples with 1 evaluation per sample. Range (min … max): 116.423 ms … 119.207 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 117.007 ms ┊ GC (median): 0.00% Time (mean ± σ): 117.182 ms ± 571.900 μs ┊ GC (mean ± σ): 0.00% ± 0.00% █ ▂ ▄▄▄▁▄▁▆▁▄▆█▆█▆▆█▁▆▆▁▄▆▁▁▁▁▁▁▆▁▄▁▁▁▁▁▁▁▁▄▁▁▄▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▄ ▁ 116 ms Histogram: frequency by time 119 ms < Memory estimate: 30.91 KiB, allocs estimate: 289.
# speedup
median(bm_cpu.times) / median(bm_gpu.times)0.7985805715942073# SP matrix multiplication on GPU
bm_gpu = @benchmark oneAPI.@sync mul!($zd, $xd, $yd)BenchmarkTools.Trial: 1 sample with 1 evaluation per sample. Single result which took 10.388 s (0.00% GC) to evaluate, with a memory estimate of 320 bytes, over 13 allocations.
For this problem size on this machine, we see GPU achieves a staggering 9 TFLOPS throughput with single precision!
# SP throughput on GPU
(2n^3) / (minimum(bm_gpu.times) / 1e9)8.467424104597974e11# SP matrix multiplication on CPU
bm_cpu = @benchmark mul!($z, $x, $y)BenchmarkTools.Trial: 1 sample with 1 evaluation per sample. Single result which took 12.125 s (0.00% GC) to evaluate, with a memory estimate of 0 bytes, over 0 allocations.
# SP throughput on CPU
(2n^3) / (minimum(bm_cpu.times) / 1e9)7.254477997545334e11We see >10x speedup by GPUs in this matrix multiplication example.
# cholesky on Gram matrix
# This one doesn't seem to work on Apple M2 chip yet
xtxd = xd'xd + I
@benchmark oneAPI.@sync cholesky($(xtxd))MethodError: no method matching UpperTriangular(::Float32) The type `UpperTriangular` exists, but no method is defined for this combination of argument types when trying to construct it. Closest candidates are: UpperTriangular(::UpperTriangular) @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:26 UpperTriangular(::AbstractMatrix) @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/triangular.jl:28 Stacktrace: [1] _chol!(A::oneArray{Float32, 2, oneAPI.oneL0.DeviceBuffer}, ::Type{UpperTriangular}) @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:188 [2] cholesky!(A::Hermitian{Float32, oneArray{Float32, 2, oneAPI.oneL0.DeviceBuffer}}, ::NoPivot; check::Bool) @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:268 [3] cholesky!(A::oneArray{Float32, 2, oneAPI.oneL0.DeviceBuffer}, ::NoPivot; check::Bool) @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:301 [4] cholesky! (repeats 2 times) @ ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:295 [inlined] [5] _cholesky @ ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:411 [inlined] [6] cholesky (repeats 2 times) @ ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:401 [inlined] [7] macro expansion @ ~/.julia/packages/oneAPI/c7eAo/src/utils.jl:71 [inlined] [8] var"##core#356"(xtxd#355::oneArray{Float32, 2, oneAPI.oneL0.DeviceBuffer}) @ Main ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:598 [9] var"##sample#357"(::Tuple{oneArray{Float32, 2, oneAPI.oneL0.DeviceBuffer}}, __params::BenchmarkTools.Parameters) @ Main ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:607 [10] _lineartrial(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters; maxevals::Int64, kwargs::@Kwargs{}) @ BenchmarkTools ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:186 [11] _lineartrial(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters) @ BenchmarkTools ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:181 [12] #invokelatest#2 @ ./essentials.jl:1055 [inlined] [13] invokelatest @ ./essentials.jl:1052 [inlined] [14] #lineartrial#46 @ ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:51 [inlined] [15] lineartrial @ ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:50 [inlined] [16] tune!(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters; progressid::Nothing, nleaves::Float64, ndone::Float64, verbose::Bool, pad::String, kwargs::@Kwargs{}) @ BenchmarkTools ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:299 [17] tune! @ ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:288 [inlined] [18] tune!(b::BenchmarkTools.Benchmark) @ BenchmarkTools ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:288 [19] top-level scope @ ~/.julia/packages/BenchmarkTools/1i1mY/src/execution.jl:461
xtx = collect(xtxd)
@benchmark LinearAlgebra.cholesky($(Symmetric(xtx)))BenchmarkTools.Trial: 1 sample with 1 evaluation per sample. Single result which took 7.649 s (0.00% GC) to evaluate, with a memory estimate of 1.00 GiB, over 3 allocations.
We don’t see GPU speedup of Cholesky at the moment.
# elementwise function on GPU arrays
fill!(yd, 1)
bm_gpu = @benchmark oneAPI.@sync $zd .= log.($yd .+ sin.($xd))
bm_gpuBenchmarkTools.Trial: 32 samples with 1 evaluation per sample. Range (min … max): 156.100 ms … 159.967 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 157.062 ms ┊ GC (median): 0.00% Time (mean ± σ): 157.240 ms ± 910.113 μs ┊ GC (mean ± σ): 0.00% ± 0.00% █ ▃ ▃█ ▃ ▃ ▃ █▁█▁▇▇▇▁██▁▁▇▁█▇▇▇▁▁█▇▁▇▁▁▇▁▁▇▇█▁▁▁▁▇▁▁▁▇▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇ ▁ 156 ms Histogram: frequency by time 160 ms < Memory estimate: 30.91 KiB, allocs estimate: 289.
# elementwise function on CPU arrays
x, y, z = collect(xd), collect(yd), collect(zd)
bm_cpu = @benchmark $z .= log.($y .+ sin.($x))
bm_cpuBenchmarkTools.Trial: 3 samples with 1 evaluation per sample. Range (min … max): 2.368 s … 2.381 s ┊ GC (min … max): 0.00% … 0.00% Time (median): 2.378 s ┊ GC (median): 0.00% Time (mean ± σ): 2.376 s ± 6.592 ms ┊ GC (mean ± σ): 0.00% ± 0.00% █ █ █ █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁█ ▁ 2.37 s Histogram: frequency by time 2.38 s < Memory estimate: 0 bytes, allocs estimate: 0.
# Speed up
median(bm_cpu.times) / median(bm_gpu.times)15.141360532518187GPU brings great speedup (>50x) to the massive evaluation of elementary math functions.
bm_cpu = @benchmark z .= tanh.($x) # on CPU
bm_cpuBenchmarkTools.Trial: 5 samples with 1 evaluation per sample. Range (min … max): 1.119 s … 1.130 s ┊ GC (min … max): 0.00% … 0.00% Time (median): 1.124 s ┊ GC (median): 0.00% Time (mean ± σ): 1.124 s ± 4.908 ms ┊ GC (mean ± σ): 0.00% ± 0.00% █ █ █ █ █ █▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁█ ▁ 1.12 s Histogram: frequency by time 1.13 s < Memory estimate: 16 bytes, allocs estimate: 1.
bm_mtl = @benchmark zd .= oneAPI.@sync tanh.($xd) # oneAPI
bm_mtlBenchmarkTools.Trial: 6 samples with 1 evaluation per sample. Range (min … max): 270.014 ms … 302.098 ms ┊ GC (min … max): 0.00% … 0.00% Time (median): 281.336 ms ┊ GC (median): 0.00% Time (mean ± σ): 283.064 ms ± 11.003 ms ┊ GC (mean ± σ): 0.00% ± 0.00% █ █ █ █ █ █ █▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁█▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁ 270 ms Histogram: frequency by time 302 ms < Memory estimate: 63.70 KiB, allocs estimate: 654.
oneAPI.jl accelerates the evaluation of tanh function by
median(bm_cpu.times) / median(bm_mtl.times)3.9936078238671047